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A RELAXATION MODEL FOR LIQUID-VAPOR PHASE 
CHANGE WITH METASTABILITY 


FRANgOIS JAMES AND HELENE MATHIS 


Abstract. We propose a model that describes phase transition including 
metastable states present in the van der Waals Equation of State. From a 
convex optimization problem on the Helmoltz free energy of a mixture, we 
deduce a dynamical system that is able to depict the mass transfer between 
two phases, for which equilibrium states are either metastable states, stable 
states or a coexistent state. The dynamical system is then used as a relaxation 
source term in an isothermal 4x4 two-phase model. We use a Finite Volume 
scheme that treats the convective part and the source term in a fractional step 
way. Numerical results illustrate the ability of the model to capture phase 
transition and metastable states. 
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!c:introduction 


1. Introduction 

In the last decades considerable research has been devoted to the simulation 
of liquid-vapor phase change, which are of major importance in several industrial 
applications. For instance liquid-vapor flows are present in water circuit of pres¬ 
surized water reactors in which the water can be submitted to saturation pressure 
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and temperature [12, 11]. The phenomena we are interested in are complex phase 
changes including the possible appearance of metastable states. An example is the 
metastable vapor which is a gaseous state where the pressure is higher than the 
saturation pressure. Such states are very unstable and a very small perturbation 
brings out a droplet of liquid inside the gas [25] . This phenomenon can appear at 
saturated pressure (or at saturated temperature for metastable liquid). It is the 
case in pressurized water reactor during a loss of coolant accident when sudden 
vaporization occurs due to the drop of pressure inside the superheated liquid [21]. 

We focus in this paper on situations where the heterogeneity of the fluid and the 
thermodynamical conditions allow the diphasic flow to be described by a compress¬ 
ible averaged model using Euler type equations. Other models can be considered, 
which account for very srnale scale by means of Korteweg type tensors including 
dispersive and dissipative effects. Such models allow to preserve metastable states 
but give only a microscopic description of the flow, see [33, 34, 1, 6, 27, 22] and 
the references cited herein. In the averaged models framework, one can distinguish 
between one-fluid and two-fluid models. 

The one-fluid model approach consists in describing the fluid flow as a single 
substance that can be present in its vapor or its liquid phase. Assuming that 
the thermodynamical equilibrium is reached instantaneously (quasi-static process), 
then the Euler system has to be closed by an Equation of State (EoS) able to depict 
either the pure phases (liquid or vapor) and the phase transition. A typical example 
is the van der Waals EoS, which is well-known to depict stable and metastable 
states below the critical temperature. However this EoS is not valid in the so- 
called spinodal zone where the pressure is a decreasing function of the density. 
This forbids the use of instantaneous kinetic exchanges, since the pressure is always 
given by the EoS, and a decreasing pressure leads to a loss of hyperbolicity in Euler 
equations, hence to instabilities and computational failure. 

To overcome this defect, and recover that the phase transition happens at con¬ 
stant pressure, temperature and chemical potential, the van der Waals pressure is 
commonly corrected by the Maxwell equal area rule [8]. This construction leads 
to a correct constant pressure in the spinodal zone but removes the metastable 
regions. 

Another way to provide a unique EoS able to depict pure phases and phase 
transition has been studied in [24, 2, 14, 17, 28]. It consists in considering that each 
phase is depicted by its own convex energy (that is its own monotone pressure law). 
According to the second principle of thermodynamics, the mixture equilibrium 
energy corresponds to the inf-convolution of the partial energies. As a result the 
mixture equilibrium energy coincides with the convex hull of the minimum of the 
two partial energies. The resulting pressure law of the mixture turns out to be 
composed of the monotone branches of the liquid and vapor pressure laws and is 
constant in the phase transition zone. Hence it is clear that such a construction 
prevents the appearance of metastable phases. 

Still in the context of one-fluid models, it is also possible to drop out the as¬ 
sumption of instantaneous equilibrium. The model involves then relaxation terms, 
which can be of various forms, but are the expression of a pulling back force to the 
equilibrium. We have to consider extended versions of the Euler system, which is 
supplemented by partial differential equations on additional quantities such as the 
volume fraction of the vapor phase or partial masses. This approach has been used 
in [5, 18, 32, 30, 26] and in [10, 15] in the isothermal case. In the later references 
the question of preservation of metastable states is adressed. 
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We mention briefly another way to describe diphasic flows, which consists in 
considering a two-fluid approach to model liquid-vapor phase change. Initially de¬ 
veloped to depict the motion of multicomponent flows [4] , such a modeling assumes 
that the fluid can locally be present under both phase. Hence the model admits 
two pressures, two velocities and two temperatures and is supplemented by addi¬ 
tional equations on the volume fraction. Phase transition is achieved by chemical 
and mechanical relaxation processes, in the limit where the kinetics is considered 
infinitely fast, see for instance [32, 35, 30]. 

One of the present drawbacks of the averaged models (one or two-fluid) with re¬ 
laxation is that there is no global agreement on the equations satisfied by the frac¬ 
tions and on the transfer terms [13]. Moreover the preservation of the metastable 
states seems to be out of reach in this framework. 

We intend in this paper to provide a model able to depict liquid-vapor phase 
change and metastable states of a single component, say some liquid in interaction 
with its own vapor. We focus on a one-fluid description of the motion while the 
phase transition is driven by transfer terms that will be coupled to fluid equations 
through a finite relaxation speed. 

The modelling of the phase transition is the core of this work. For the sake 
of simplicity we assume the system to be isothermal. We propose transfer terms 
obtained through the minimization of the Helmholtz free energy of the two-phase 
system. We use for both phases the same equation of state which has to be non 
monotone, typically the reduced van der Waals equation. This choice allows us to 
recover all possible equilibria: pure phases (liquid or vapor), metastable states and 
coexistence states characterized by the equality of pressures and chemical potentials. 
These are physically admissible states, but the set of equilibrium states also contains 
the spinodal zone, which is irrelevant. Thus the key point now is to characterize 
the physical stability, and hence admissibility, of these states. It turns out that this 
has to be done in terms of their dynamical behaviour. More precisely, we design a 
dynamical system which is able to depict all the stable equilibria of the system as 
attraction points. In particular we show that metastable states and mixtures have 
different basins of attraction, so that they can be differentiated only by their long¬ 
time behaviour with respect to initial conditions. Hence there is no hope to recover 
both metastable states and coexistent states under the assumption of instantaneous 
equilibrium, which amounts precisely to choose a priori this long-time behaviour. 

This dynamical system is used as a transfer term in an isothermal two-phase 
model in the spirit of [29] and [3]. The extended Euler system we obtain in this 
way provides a regularization of the isothermal Euler equation with van der Waals 
EoS, which takes the form of a mixture zone surrounding the physical interfaces, 
see Section 5 below. 

The paper is organized as follows. Section 2 is devoted to the thermodynamics 
of a two-phase fluid. We provide the definitions of the common potentials and give 
some details on the reduced van der Waals model. Assuming that both phases 
follow the same non monotone EoS, we describe the thermodynamical equilibrium 
as the result of a minimization process on the Helmoltz free energy of the two- 
phase fluid. The section ends with the study of the equilibria of the optimality 
system. Section 3 is the core of this work. It is devoted to the construction of the 
dynamical system based on the results of the previous section. A few numerical 
simulations illustrate the ability of the system to catch both the Maxwell line and 
the metastable states in the van der Waals EoS. The dynamical system is then 
plugged as a source term in a 4 x 4 isothermal model in section 4. We provide a 
study of the homogeneous system, which is conditionally hyperbolic. However we 
prove that for smooth solutions the hyperbolicity regions are invariant domains of 
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the system with relaxation. In section 5 we present several numerical illustrations 
which assess the ability of the model to deal with metastable states. They are 
obtained using a classical finite volume schemes that treats the convective part and 
the source terms with a time-splitting technique. 


2. Thermodynamics and the van der Waals EoS 

2.1. Thermodynamics of a single fluid. Consider a fluid of mass M > 0 occu¬ 
pying a volume V > 0, assumed to be at constant temperature T. If the fluid is 
homogeneous and at rest, its behavior is entirely described by its mass, its volume 
and its Helmholtz free energy F. According to Gibbs’ formalism [16], the fluid is 
at equilibrium if its Helmholtz free energy is a function, also denoted by F, of its 
mass M and volume V: 

(1) F : (M, V) —> F(M, V). 

Notice that we do not address yet the stability of equilibrium states. At this 
level the involved quantities are extensive, which means that they share the same 
scalings as the volume V. This corresponds to the notion of homogeneous sample: 
for twice the volume, the mass is double, and the energy as well. The mathematical 
consequence of this notion is that extensive variables have to be related through 
positively homogeneous functions of degree 1 (PHI), namely 

(2) VA > 0, F(XM, XV) = XF(M, V). 

We assume in addition, and without loss of generality, that the energy function F 
belongs to C 2 (R + x R + ). 

Remark 2.1. The regularity of F seems to preclude phase transitions, but this is 
not the case because no convexity assumption is made at this stage. We shall come 
back to this in more details in the next section. 


We introduce now intensive quantities, by which we mean that they do not 
depend on the volume scaling. A typical example is the density p = M/V, but such 
quantities also appear as derivatives of the the equilibrium relation (2), which are 
homogeneous of degree 0 by construction. In this way two fundamental quantities 
are to be considered, namely the pressure p and the chemical potential p, defined 

by 


( 3 ) 


dF 

P = -^(M, n 


dF 


Notice that these quantities are defined only when the system is at equilibrium, 
and we recover the classical thermodynamic relation for isothermal flows 


( 4 ) 


dF = —pdV + pdM. 


Another classical property of thermodynamical potentials is the so-called Gibbs 
relation, which results from the Euler relation for PHI functions: 

F(M,V) = VF(M,V)-(*f) . 


Using (3), we obtain 

(5) F(M, V) = pM - pV. 


For most of the following computations it is useful to rewrite the preceding 
relations using intensive variables. For a fixed volume V, we denote / the Helmoltz 
free energy per unit volume, which is a function of the density p = M/V: 
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( 6 ) 


f(p) = f 





We keep the notations p and p the pressure and the chemical potential as functions 
of the density p: 


(7) p(p)= p (^pj 


dE (M \ 

W [v' 1 ) ’ 


p{p) = p 



dE /M \ 
dM V^’7 ' 


Thus we obtain an intensive form of the Gibbs relation (5) 


( 8 ) 


f{p) = pp{p)~p{p)- 


On the other hand, from the definitions of p, p and / we obtain easily the following 
relations 


(9) p{p) = f'(p ), p(p) = pf(p) - f(p), pp'(p) = p\p). 


2.2. Example: the van der Waals EOS. The extensive energy of a van der 
Waals (monoatomic) fluid is 

(10) F[M , V) = -^+Rt(m log ^j^-M^, 


where R stands for the perfect gases constant and a and b are positive constants 
(a = b = 0 leads to the perfect gases law). Since we consider the isothermal model, 
the temperature T is a parameter here. We refer to [25, Ch. 7] for a justification 
of this law from statistical thermodynamics. The constant b is proportional to the 
proper volume of a particule such that V > Mb, and the potential — depicts the 

interaction between particules. 

The pressure and the chemical potential associated read 


( 11 ) 



| M 
V 2 V-Mb ’ 



+ RT log 


M 

V-Mb 


+ RT 


Mb 

V-Mb' 


The intensive quantities are 


f(p) = -ap 2 + RTp (log - l) , 

(12) p(p) = ~ap 2 + { >m i)p , 

,.(P) = -up + RT log + ■ 

The behavior of the isotherm curves in the pressure-density plan is depicted in 
Figure 1. The critical temperature T c is the lower limit of temperatures such that 
the pressure is an increasing function of the density. At T = T c the pressure curve 
admits an unique horizontal inflection point, called the critical point, denoted C 
on Figure 1. 

In the sequel we will consider the reduced form of the van der Waals equation 
(see [25, Ch. 8]). Denoting ( p c ,p c ) the coordinates of the critical point one has 


p'{pc) = -‘ZaPc + 


RT C 


(1 - bp c ) 2 


= 0, 


"/ \ o i 2 bRT c 

p M = - 2a + jrr^j 3 =°- 
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fig:vdw_isoth 


Figure 1 . Phase diagram for the van der Waals EoS in the (p , p) 
plan. Below the critical temperature Tq, the isotherm curve de¬ 
creases in the spinodal zone which is delimited by the densities 
p~ < p + . In that area the isotherm is commonly replaced by an 
horizontal segment (dashed line) which coincides with the isobaric 
line at constant pressure p = p*. Such a construction defines the 
two densities p\ and p 



Considering normalized critical quantities, that is setting p c = 1, p c = 1 and T c = 1, 
one obtains the reduced van der Waals EOS with R = 8/3, a = 3, and 6 = 1/3 and 


eq:reduced_vdw 


f(p) = -3 P 2 + |P (log(3p/(3 - p)) - 1), 

(13) p(p) = -3p 2 + 

o p 

P(p) = -6 P + log(3p/(3 - p)) + |(3 p /(3 - p)) 

for 0 < p < 3. 

Below the critical temperature T c the pressure is not monotone with respect 
to the density (see the red curve on Figure 1): in a region called spinodal zone, 
delimited for a given temperature by the densities p~ < p + , the pressure decreases 
with respect to the density and thus leading to unstable states. According to (9) 
the energy / is non-convex in the spinodal zone. 

In that region the isotherm is commonly replaced by the Maxwell area rule in 
order to recover that phase transition happens at constant pressure and constant 
chemical potential. The Maxwell construction is commonly applied on the pressure 
(see [8] for instance) in such a way that the two zones delimited by the van der Waals 
isotherm and the Maxwell line (above and below the Maxwell line respectively) have 
the same area. This is not the case in Figure 1 because in our context, the equal 
area rule is obtained on the chemical potential, see Section 2.4 below (Proposition 
2.1). In any case, the idea is that the isotherm curve is replaced locally by the 
horizontal segment, the so-called Maxwell line, that coincides with some isobaric 
line p = p*. Such a construction defines two densities, denoted p\ and p\, as well 
as the value of p*. An equivalent way to compute this correction is to build the 
convex hull of the function /, see [2] and [17]. 

However this construction removes the admissible regions delimited by the densi¬ 
ties where the pressure law is still nondecreasing. Such regions are called metastable 
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regions. At a given temperature these regions correspond to densities belonging to 
the range [p\, p~] and \p + ,p 2 }. 


2.3. Thermodynamics of phase transition. The van der Waals model depicts 
the thermodynamical behavior of a pure substance under its liquid state, gaseous 
state and the coexistence state. The non convexity of the EoS allows to capture 
metastable states but does not give a relevant representation of the coexistence 
phase. A convenient way to cope with this problem is to represent the fluid under 
consideration as a set of several copies of itself under different pure phases (liquid 
or gaseous phases). Such a representation is used in [24, 5, 18, 2, 32, 17, 14, 30, 26] 
assuming that the fluid is described by two copies, each one satisfying a convex EoS 
that differs from the one of another copy. See also [23] where such a representation 
is used in the context of adsorption-desorption of a mixture. 

We adopt here a slightly modified approach. Let us consider I copies of the 
pure substance, I > 1 being some integer a priori arbitrary. Each copy is depicted 
by its mass Mi > 0 and its volume V) > 0, and is assumed to be at equilibrium. 
Hence a copy can occupy a volume with zero mass. We suppose that each copy 
follows the same non-convex energy function F(M l: V)), typically the van der Waals 
extensive energy given by (10). This assumption contrasts with the aforementioned 
references where different convex energy functions are considered, and the number 
of copies is prescribed. 

Thanks to the mass conservation, the complete system has a total mass M = 
Xa=i Afj. Assuming that the copies are immiscible, the total volume is V = Xu=i 
(for a mixture of gas, one has V) = V, Vi = 1,... I, and this condition implies 
Dalton’s law, see [8]). Out of thermodynamical equilibrium the free Helmoltz energy 
of the system reads 

i 

=5Z*W,Vi). 

i—1 

Let us fix the total mass M and volume V of the system. According to the 
second principle of thermodynamics (see [8] for instance), the thermodynamical 
equilibrium corresponds to the solutions of the constrained optimization problem 


inf 

I>l,Mi>0,Vi>0 


{^2 F(M u Vi)-, ^ Vi 


I 

v, J2 M > = M J- 

i— 1 


We stress the fact that the total number I of possible copies is not fixed a priori 
here. However, as a consequence of Caratheodory’s theorem in dimension 2 (see 
e.g. [31, Ch. 17]), we can state the Gibbs phase rule (see [25, Ch. 9]), which gives 
the expected result in this context. 


Lemma 2.1 (Gibbs phase rule). We have I < 2. 


In the sequel we use just as well the term phase as the term copy. 

Taking into account Lemma 2.1 and the above intensive formulations, we rewrite 
the constrained optimization problem for a fixed mass M and volume V, hence for 
a fixed p, as 

(14) inf Jaif(pi) + a 2 f(p 2 )}, 

P i>0 

p 2 >0 

under the constraints 


(15) a\+a 2 = 1, 

(16) aipi+a 2 p 2 = p. 
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Here an = -p £ [0,1] denotes the volume fraction and pi = ^ > 0 the density of 

the phase i = 1,2. Notice that we removed the optimization on the phase number 
I, since single phases can be recovered by p\ = P 2 (and any af) or one a, = 0 (with 
undetermined p : j, j f i). 

We rule out the equality case by noticing that, provided pi 7 ^ p 2l we can rewrite 
the constraints (15)-(16) as 

xR + -> [0,1] a 2 :l + xR + xI + -> [0,1] 


ai : R+ x 


(17) 


(p,Pl,P2) 


P2 


Pi - P2 


(p,Pl,P2) -- 


Pi 


Pi - P2 


We have or > 0 if and only if pi < p < P 2 or pi > p > p 2 . Therefore, accounting 
on the reduced form of the van der Waals model (13), we shall assume in the sequel 
that the densities p, p\ and P 2 satisfy 


(Hi) 


0 <pi < p < p 2 < 3 and Pi < p 2 - 


This is not a restriction, as we shall see below (see Proposition 3.1). 

Mi 

One can also introduce the mass fraction ipi = i = 1,2 defined by 


Vi : 


(18) 


(P,Pl,P2) ^ -J 


+ ^[ 0 , 1 ] 

1 1 

P2 


T2 : 


[ 0 , 1 ] 

1 


1 


(P,Pl,P2) ^ -J- 


1 


1 

Pi 

1 ’ 


Pi P2 P2 Pi 

that satisfy P 1 + P 2 = 1 and <pi L P 2 > 0 if and only if the assumption (Hi) is satisfied. 
Such quantities will be useful in the mathematical study of the isothermal two-phase 
flow model introduced in Section 4. 


eq:infconv 


Remark 2 . 2 . In the aforementioned references [24, 5, 18, 2, 17, 14] the method 
consists in describing the two-phase fluid by a coexistence of two copies of the same 
substance. The description can be made either on the extensive variables or the 
intensive one. Unlike our present approach the two copies do not follow the same 
EoS: each copy is described by its own strictly convex extensive energy Ft, i = 1,2, 
which is a function of the mass Mi and the volume Vi of the phase. Following the 
second principle of thermodynamics [8] , at equilibrium, the extensive energy of the 
two-phase fluid is given by 

(19) F((M, V)) := FiUF 2 (M 1 V) = min Fi(M u Vi) + F 2 (M - M u V - Hi), 


under the constraints of mass conservation M = Mi + M 2 and immiscibility (with¬ 
out vacuum) V = Vi + V 2 . This operation □ is called inf-convolution operation in 
the convex analysis framework [19]. In [17] the authors investigate the link between 
the inf-convolution, the Legendre transform and the (max, +) algebra. The Legendre 
transform of a energy F is a convex function (M*,V*) —X F*(M*,V*) defined by 

F*(M*,V*)= inf {M*M+ V*V - F(M,V)}. 

M*> 0,V*>0 


The inf-convolution is transformed into an addition by the Legendre transform which 
implies that 

(FiDF 2 )* = F* + F*. 

Moreover in the case of convex lower semi continuous (sic) functions Fi, i = 1,2, 
one has 


FiDF 2 = (FiDF 2 )** = (F* + F 2 *)*. 
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It means that the energy F eq of the two-phase fluid at equilibrium is given by 
(20) F eq = F **, 

where F** is the convex hull of the energy F. As it is proved in [2] the construction 
of the convex hull of the energy F is equivalent to the Maxwell construction. Hence 
the operation ( 20 ) removes the metastable regions. 

2.4. Equilibrium states. This section is devoted to the characterization of the 
equilibrium states of the thermodynamical system, that is states that realize the 
inhmum in the optimization problem (14)-(15)-(16). The function we minimize is 
defined by: 

f : [0,1] x [0, l]xK + xl+-tl 

(ai, a 2 , pi, p 2 ) i-t «i/(pi) + a 2 f{p 2 ), 

and is C 1 on the space {pi > 0, i = 1,2}. The constraints (15)-(16) are affine, 
so they are also C 1 . We are thus in position to use the Lagrange multipliers 
characterization of the infimum (see [31, Ch. 28]): X a £ K. and X p £ R. respectively 
correspond to the constraints (15) and (16). 

Using the definition ( 8 ) of the free energy / and the pressure and chemical 
potential definitions (9), one deduces the following optimality system of equations 


( 21 ) 

f{pi) + Aq, + XpPi 

= 0 , 

( 22 ) 

f{P2) + Aq, + \pP 2 

= 0 , 

(23) 

aip{pi) + Apoq 

= 0 , 

(24) 

a 2 p{p 2 ) + X p a 2 

= 0. 


From this optimality system we deduce immediately that there are two different 
kinds of equilibria. 

Lemma 2.2. Under hypothesis {Hi), the equilibrium states are 

(1) Pure liquid or gaseous states: or = 0 (resp. a 2 = 0), with p 2 = p, 
pi < p arbitrary (resp. pi = p, p < p 2 arbitrary) 

(2) Coexistence states: aia 2 ^ 0, with (pi,p 2 ) satisfying p{pi) = p{p 2 ) and 
P(Pl) =P{P2)- 


Proof. The case or = 0 corresponds to the saturation of the constraint or + a 2 = 1 
see (15). It leads to a 2 = 1 and thus p- 2 = p that is only the phase 2 is present. 
Conversely if a 2 = 1 only the phase 1 remains. 

On the other hand let us assume aia 2 ^ 0. Then (23) and (24) lead to the 
equality of the chemical potentials 

P{Pl) = P{P2) = -A p . 

Then the intensive Gibbs relation (8) allows to rewrite the conditions (21) and (22) 
as 

PiP{Pi) - p(Pi) + K + A P pi = 0, * = 1,2. 

Since —A p = p{pf), this leads to the pressures equality 

P(Pl) =P{P2) = A a . 


□ 


To proceed further we need the following result for coexistence states. 

Proposition 2 . 1 . Under hypothesis {Hi) and if aia 2 7 ^ 0, the following proposi¬ 
tions are equivalent and uniquely define the couple of densities p\ < p 2 . 
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(1) The chemical potentials and the pressures are equal 

(25) p(p\) = p{pV>, 

(26) P{P*)=P{P*2)- 

(2) The Maxwell’s area rule on the chemical potential holds 

(27) [ n{p 2 +t(p! - p 2 ))dt = p{p\) = p{p* 2 ). 

Jo 

(3) The difference of the Helmoltz free energies reads 

(28) f{p 2 ) - /(Pi) = p(pl)(p*2 ~ P*i) = p{p* 2 ){p*2 - Pi)- 

Proof. The identities (27) and (28) are equivalent since f'{p) = p(p), see (9). Now 
assume (25)-(26) hold. Then the intensive Gibbs relation (8) gives (28). Conversely 

(28) implies the chemical potentials equality and thus the pressures equality. Now 

the uniqueness of (Pi,p 2 ) follows easily form a geometrical argument using the 
Maxwell’s area rule. □ 

The most famous characterizations of diphasic equilibria are (25)-(25) and (27), 
although the latter is usually written in terms of pressure. We can recover this 
form by writing the intensive relations with the specific volume t = V/M instead 
of p. The third relation (28) is not so classic but will be useful in the sequel. 

The density p\ (resp. p 2 ) separates the range of pure gaseous state and the range 
of metastable gas (resp. the range of pure liquid state and the range of metastable 
liquid), see Figure 1. These densities define the coexistence pressure p* and the 
coexistence chemical potential p*: 

(29) p* = p(pl) = p{p* 2 ), p* = p(pl) = p(p 2 ). 

We emphasize that the necessary conditions for equilibrium contain both unsta¬ 
ble (spinodal zone), metastable and stable states. Nothing at this stage can make 
the difference, which turns out to be of dynamical nature, in a way we make precise 
now. 


3. Dynamical system and phase transition 

This section is devoted to the construction and the analysis of a dynamical 
system deduced from the optimality conditions given in Lemma 2.2 and Propo¬ 
sition 2.1. We will prove that the equilibria of this dynamical system are either 
pure liquid/vapor states, pure metastable states or coexistence states in the spin¬ 
odal zone (that is states satisfying the properties (25)-(27) of Proposition 2.1). We 
emphasize that the difference between metastable states and coexistent states actu¬ 
ally relies on the long-time dynamical behaviour of the solutions to the dynamical 
system. No static characterization can be given. 

The section ends with numerical illustrations that assess the ability of the dy¬ 
namical system to preserve metastable states. 

3.1. Construction of the dynamical system. We want to construct a dynami¬ 
cal system which equilibria coincide with the equilibria of the optimization problem 
depicted in Lemma 2.2. To do so, we force the dynamical system to dissipate the 
Helmoltz free energy defined by 

J:M + xR + xK + aR 

(p,pi,P 2 ) a 1 (p,p 1 ,p 2 )f{pi) + a 2 (p,pi,P 2 )f{p 2 ), 

under the optimization constrains (15) and (16). 

Assuming that p, p\ and p 2 are now time-dependent functions, we can compute 
the derivative of the total Helmoltz free energy T with respect to time and deduce 
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eq: r 


eq:gradFl 


eq:gradF2 


eq:relat_f 


eq:Fdot 


eq: syst_dyn2 


appropriate time derivatives of p , p\ and p 2 such that T is dissipated in time. 
Moreover we want pure states (either liquid, vapor or metastable) to be equilibria 
of the dynamical system. Hence we have forced the time derivatives of p, pi and 
P2 to vanish in case of pure state (that is when ol\Ol 2 = 0). 

For sake of readability we denote 



the vector of admissible densities satisfying the assumption {Hi). By the defini¬ 
tion (17) of the volume fractions a;, i = 1,2, one has 


V r ai(r) = — V r a 2 (r) = 


From this and the definition (30) of J -, we easily get the gradient of the free energy 



(32) V r J 7 (r) = 


f{pi) - f(P2) 

ai{r){p 2 {p{p 2 ) ~ p(pi)) +p(pi) ~P(P 2 )) 
K a 2 {r){pi{p(p 2 ) - p{pi)) +p{pi)~p{p 2 ))-j 

Note at once that it can be expressed in terms of relative free energy 


1 

Pi - P2 


(33) 


Vr^(r) = 


1 


( f(pl) - f(P2) 
a 1 {r)f{p 2 \pi) 


Pl P2 \-a 2 (r)f{p 1 \p 2 )j 

where the relative free energy of p 2 with respect to p\ is defined by 


(34) f {p 2 \pi ) := f{p2) - f{pi) - p{pi)(P2 - Pi)- 


We turn now to the definition of our dynamical system. Because of the mass 
conservation, we obviously impose that 


P = o. 

The main idea to proceed further is that we want the system to dissipate the total 
Helmoltz free energy T. Combining the definition (30) of F and the expression of 
V r J r (r), one computes the time derivative of F along some trajectory, that is 

(35) F{r) = — 1 — [a\{v)f {p 2 \pi)pi - a 2 {v)f{pi\p 2 )p 2 ). 

Pl - P2 

We now propose the following dynamical system 

[p= o, 

(36) < Pi = {p Pi){p P 2 )f {p 2 \pi)i 

[P2= +{P- Pl){P~ P2)f{Pl\p 2 )- 

An easy computation shows that this system dissipates J- along its trajectories (un¬ 
der the assumption {Hi)). Indeed using the expression (35) of T and the equations 
of pi and p 2 , one gets 

•^(r) = ~{P - Pi) [ a i{ r )f{P 2 \pi)] 2 + {p~ Pi) [a 2 {v)f{pi\p 2 )f . 

The multiplicative term {p— pi){p — p 2 ) in the equations on pi and p 2 ensures that 
the right hand side of the system vanishes in case of pure states (either pure liquid, 
vapor or metastable states). 


Remark 3.1. We emphasize that the choice of the right hand side of (36) is 
somewhat arbitrary. Other terms might be more efficient, but this one was chosen 
for its simplicity and its interpretation in terms of relative free energy. 












eq:dyn_syst_equ 


sec:study-dynam-syst 


prop:DynProp 


eq:Hl_time 


eq:alphal_time 


th:attraction 
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One can check easily that a\p\ + a 2 p 2 = 0 so that it is consistent with the total 
mass conservation equation p = 0 . 

An equivalent dynamical system can be written in terms of the time derivatives of 
the volume fractions and of the partial masses aipi, i = 1,2. Accounting for the 
constraints (15)-(16) and for the system (36), some straightforward computations 
lead to 


(37) 


Oil 


Oil Pi 


= —0:2 

= Ot\(p - Pl)f(p 2 \pi) + Oi\{p - p 2 )f(pi\p 2 ), 

= — 02^2 

= Oi\p 2 {p - pi)f(p 2 \pi) + a\pi(p - p 2 )f(pi\p 2 ). 


This formulation does not allow to compute single phase systems, since the deter¬ 
mination of the partial densities pi is impossible when o^ = 0. Thus we rather use 
the dynamical system (36). 


3.2. Equilibria of the dynamical system. The major result of this paragraph 
is that the attractors of the system are either pure liquid/vapor states, including 
metastable states, or the coexistence state defined by (27)-(26), see Proposition 2.1. 

Proposition 3.1. The dynamical system (36) satisfies the following properties. 

(1) If the assumption (Hi) is satisfied at t = 0, then it is preserved in time (in 
particular, the assumption p± < p 2 is preserved). 

(2) If oi(0) = 0 (resp. Oi(0) = 1) then oi (t) = 0 (resp a±(t) = 1), for all 
time. 

Proof. First we address the preservation of the hypothesis (Hi). Some straightfor¬ 
ward computations lead to 

(38) pi - p 2 = (p - pi)(p - p 2 )(pi - p 2 )(p(p 2 ) - p(pi))- 

Thus if pi — P 2 is zero at t = 0 then this property is preserved for all time and the 
sign of the difference pi — p 2 is also preserved. The property on the volume fraction 
a.i is proved in the same way noting that 

(39) di(r) = —ai(r)(p - pi)(pi(p(pi) - p(p 2 )) - p(pi) +p(p 2 ))- 

□ 


We turn now to the characterization of the equilibria of the dynamical system. 
Since the total mass p is constant in time, they are parametrized by p. Moreover it 
arises that their thermodynamical characterization can be given in terms of attrac¬ 
tion basins. We refer the reader for graphical references to figures 2, 3 and 4, where 
attraction basins are drawn as hatched zones. The last two figures correspond to 
the gaseous phase, similar pictures can be drawn for the liquid phase. 

Theorem 3.1. Assume that the initial data (p(0),pi(0),p 2 (0)) of the system (36) 
satisfy (Hi). Then the equilibria are given by 

(1) spinodal zone: p\ < p( 0) < p\. The unique equilibrium is ( p\,p 2 ) char¬ 
acterized by Proposition 2.1, with ai given by (17). The attraction basin is 
(0,p) x (p, 3). 

(2) pure gaseous states: p(0 ) < p\ (resp. pure liquid states: p(0) > pt,). 
The set of equilibria is {p} x (p, 3), with 0 ^ = 1 (resp. (0,p) x {p} with 
a 2 = l). The attraction basin is (0 , p) x (p, 3). 

(3) metastable states: p\ < p( 0) < p~ (resp. p + < p(0) < p\). There are 
two sets of equilibria 
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(a) perturbation within the phase: ^ 2 ( 0 ) < p~ (resp. pi(0) > p + ). 
The set of equilibria is {p} x (p,p~), with or = 1 (resp. (p + , p) x 
{p} with Oii — 1/ The attraction basin is then (0 , p) x (p,p~) (resp. 
( P + ,P ) x (P, 3)j- 

(b) outside perturbation: p~ ' < P ‘1 (0) < P* 2 (resp. p + > pi(0) > 
There is a unique equilibrium (p\,p 2 ), characterized by Proposition 
2.1, with oii given by (17). The attraction basin is (0 ,p) x (p~, 3) 
(resp. (0,p+) x (p,3)j. 


Proof. We look for Lyapunov functions for each case. 
Spinodal zone. Let us define 


eq:SpinLyap 


eq:PureLyap 


(40) G(r) = aif(pi) + a 2 f(p 2 ) - p*(cnpi + a 2 p 2 ) + p*(ai + a 2 ), 

where p* and p* are defined by the Maxwell area rule (29) on the chemical potential 
p, and ai are the functions of p,p\,p 2 given by (17). The intensive version of the 
Gibbs relation ( 8 ) implies that G(r*) = 0 where r* = (p, p\, p 2 ) T ■ Straightforward 
computations show that V r G = V r T — (p*,0, 0) T , this implies that V r G(r*) = 0, 
and that G(r(f)) < 0 since the free energy is dissipated. Hence the function G is a 
Lyapunov function on the admissible domain defined by (Hi), which contains the 
unique equilibrium r*, see Figure 2. 

Pure stable states The expected equilibrium states are now r = (p,p,p 2 ) for 
any Pi > P, see Figure 3. We introduce the function 

(41) G ( 1 ) (r) = f(pi) - ppi +p, where p. = p(p), p = p(p). 

Once again, the Gibbs relation leads to G^^r) = 0 for any equilibrium r, and we 
easily have V r GG)(f) = 0. Now we have 


eq:IneqConv 


^G (1) (r(f)) = -(p(p-i) - p(p))(p - pi)(p- p 2 )f(p 2 \p 1 ). 

Since p\ < p < p 2 and p is nondecreasing on ]0,pj] the right-hand side in the 
preceding relation is nonpositive if /(P 2 IP 1 ) > 0. But using again the Gibbs relation 
and (34), this can be rewritten 

(42) f(pi) - f(pi) > n(pi)(pi -pi). 

This convexity inequality holds true for all (pi,p 2 ) such that p\ < p\, by the very 
definition of p\ (see the definition (20) of F**). Since we consider pure liquid states, 
one has p\ < p < p\. Hence GW is dissipated in time, leading to the conclusion. 

Metastable states Two types of equilibria are encountered in this situation, 
with two distinct attraction basins, see Figure 4. 

• the metastable basin, which appears with \\\ hatches in the figure, corre¬ 
sponds to perturbation of a given state within the same phase. It is actually 
an attraction basin, because the function G^ 1 ^ defined above (41) is also a 
Lyapunov function in this domain. Indeed the convexity argument still 
holds true for any (pi,p 2 ) such that p 2 < p~, since / is convex below p~. 

• the unstable basin, corresponding to perturbations of the state by the other 
phase (III hatched zone). This basin is governed by the same Lyapunov 
function as for the spinodal zone (40). 

□ 


Remark 3.2. This is a formalization of the remark in Landau & Lifshitz [25, §21] 
“[...]we must distinguish between metastable and stable equilibrium states. A 
body in a metastable state may not return to it after a sufficient deviation” 
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fig:SpinodalBasin 



Pi 


Figure 2. Spinodal states. The blue area refers to nonattainable 
states according to {Hi). The attraction basin is the hatched area 
(///). The unique attraction point (p*, p 2 ) appears in red. The 
green zone corresponds to the invariant hyperbolicity region (see 
Section 4.4). 


fig:PureBasin 



Figure 3. Pure gaseous states. The blue area refers to nonattain¬ 
able states according to (H 1 ). The attraction basin is the hatched 
area (\\\). The set of attraction points appears in red. The two 
green zones correspond to the invariant hyperbolicity regions (see 
Section 4.4). 


sec: resol-immer 

3.3. Numerical illustrations of the dynamical system behavior. We pro¬ 
vide in this paragraph some numerical examples to illustrate the behavior of the 
dynamical system (36). First we give the characteristic values of the reduced van 
der Waals EoS at a fixed nondimensionalized temperature T = 0.85. Figure 5 
presents the corresponding isothermal curve in the (p, p) plan. 

The densities p~ and p + correspond to the extrema of the chemical potential. 
The densities p\ et p 2 are obtained by the Maxwell’s equal area rule construction 
on the chemical potential (27) or equivalently by solving (25)-(26). 
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fig:MetaBasin 



Figure 4. Pure gaseous metastable states. The blue area refers 
to nonattainable states according to (if i). The attraction points 
appear in red. The corresponding attraction basins are the hatched 
areas. For any state belonging to the hatched area (\\\), the at¬ 
traction points are metastable vapor sates such that {pi = p}. For 
any state belonging to the hatched area (///), the attraction point 
is the coexistence state The two green zones correspond 

to the invariant hyperbolicity regions (see Section 4.4). 


Chemical potential - density plane: isothermal curves 



fig:mu_T=0.85 


Figure 5. Isotherm curve in the (p, p ) plan at T = 0.85. 


Table 1 contains the values of the characteristic densities and the corresponding 
values of the pressure and chemical potential at T = 0.85. 

Using a Backward Differentiation Formula (provided by Scilab for stiff problem), 
we provide numerical illustrations of the attraction basins depicted in Theorem 3.1. 
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p 

P(P ) 

p(p) 

p- = 0.5810799446067 

0.62055388470356498 

-3.68447967881140137 

p+ = 1.4888047089018 

0.04962960899844759 

-4.24339302065563029 

pi = 0.31972996451885 

0.504491649787487 

-3.97717851100986 

p* 2 = 1.8071403273364 


Table 1. Characteristic densities at T = 0.85 and corresponding 
pressure and chemical potential values 


tab:carac 


To do so, we fix 50 initial values of p(0) in a given interval and set the partial densi¬ 
ties pi( 0) as perturbations of p(0). Then the dynamical system is solved for a final 
time Tf = 10 3 with a time step set to 10 _3 s. We provide the graphs of the quanti¬ 
ties a\ (Tf), Pi{Tf) and a\(Tf)p(p\{^Tf)) + ct 2 (Tf)p(p 2 (Tf)), plotted with respect to 
the density p(0), for each one of the 50 initial conditions (p(0), pi(0), P2(0)). The 
mixture pressure profile is compared with the classical van der Waals pressure curve 
and the Maxwell line. 

Metastable states: The initial density p(0) takes on 50 values in [p\, p~], For 
this computation, we set pi(0) = p(0) — 0.1 and P 2 ( 0 ) = p(0) + 0.2 so that we can 
observe the perturbation within and outside the metastable vapor zone. According 
to Figure 6-top left, the mixture pressure presents two different parts: the left part 
(for p < 0.45) remains on the van der Waals pressure curve and a second part 
coincides with the Maxwell line. The first part corresponds to the perturbation of 
the metastable vapor state within the phase, while the right part corresponds to 
the perturbation outside the metastable vapor, that is when p~ < p2(0) < p*%- One 
can check on the volume fraction curve (see Figure 6-top right) that a\(Tf) = 1 
for p(0) < 0.45, that is only the phase 1 is present. Then 0 < or < 1 and the 
corresponding final partial densities Pi(Tf ), i = 1,2, coincide with the densities p\ 
and P2 respectively, which explains that the pressure matches with the Maxwell 
line. 

Perturbation of the density p in the whole domain: Figure 7 corresponds 
to an initial density p( 0) which takes on 50 values between 0.2 and 1.8, while 
pi(0) = p(0) —0.1 and ^ 2 ( 0 ) = /o(0) + 0.1. One observes that an initial perturbation 
of the density p(0) leads to final states which belong to either pure vapor/liquid 
states, including metastable states, or the coexistent state. Hence the mixture 
pressure coincides with the admissible branches of the van der Waals pressure curve 
or with the Maxwell line. 

However the convergence is not obvious for p(0) close to p~ on the left (resp. to 
p + on the right). This can be observed in Figure 7 top-right (plot of the volume 
fraction), where the parts on the left and on the right should be straight lines. 
Actually, since for p close to p~ the perturbation chosen is in the spinodal zone, 
we expect the equilibrium to be on the Maxwell line, which is not the case. We 
suppose that the final time is not large enough to ensure the actual convergence. 
As a matter of a fact, we observed that the requested time to reach convergence is 
larger in the metastable zone than in the spinodal one. 


sec:model 


4. The isothermal model 

This section is devoted to the definition and study of a 4 x 4 van der Waals 
isothermal two-phase flow model. Since we are interested in the modeling of phase 
transitions with possible metastable states, the liquid-vapor flows that we consider 
are submitted to strong thermodynamical perturbations. Hence we propose to 
depict the dynamic of the flows by a compressible averaged model, namely Euler 
type equations. 
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fig:sci_pert_meta 


Figure 6. Numerical illustration of a perturbation within the 
metastable vapor zone. The initial density p( 0) takes on 50 val¬ 
ues in [p\,p~] and Pi(0) are perturbations of p( 0) under the as¬ 
sumption (Hi). The top-left figure corresponds to the mixture 
pressure at final time Tf = 10 3 . For densities p(0) < 0.45, the 
pressure coincides with the reduced van der Waals pressure, while 
for p(0) > 0.45 it matches with the Maxwell line. One notices 
that the volume fraction a±(Tf) is either constant equal to 1 (for 
p(0) < 0.45), correspondong to pure phase 1, or takes on values in 
]0,1[, which means that the system reached the coexistent state. 


In order to model phase transitions, the hydrodynamic part of the model is 
classically coupled with a relaxation source term which carries on the mass transfer. 
Since we wish to take into account possible metastable states, the equilibria of the 
source term have to be either pure liquid/vapor states, metastable states or the 
coexistence state given by (27)-(26). Hence we propose a coupling between the 
dynamical system (36) introduced in the previous section and a modified version 
of the isothermal two-phase model proposed in [3]. 

After defining the model, we study several properties of the system, such as ex¬ 
istence of a decreasing energy, hyperbolicity and Riemann invariants for the homo¬ 
geneous system. Notice at once that we have only partial results for hyperbolicity, 
as noticed before in the literature, because of the spinodal zone. This leads to a 
formal study of invariant hyperbolicity domains in the last subsection. 

4.1. Definition of the model. The basic isothermal Euler system contains the 
balance equations accounting for the conservation of the total mass and the total 
momentum of the two-phase flow. We propose to extend this system with two 
equations describing the evolution of the partial densities p\ and p 2 which are now 
functions of time t and space x. The two phases evolve with the same velocity u. 
The momentum equation involves a pressure term which is the mixture pressure 
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fig:sci_meta 


Figure 7. Numerical illustration of a perturbation of p(0) in 
the whole domain. The mixture pressure a\(Tf)p(p\(Tf)) + 
a 2 (Tf)p(p 2 (Tf)) coincides with the admissible branches of the re¬ 
duced van der Waals pressure in the pure liquid/vapor states, in¬ 
cluding metastable state. In the spinodal zone it corresponds to 
the Maxwell line. 


a ip(pi) + cx 2 p(p 2 ). Here a\ and «2 are given by (17) but for the sake of readability 
we skip this dependence in what follows. 

The system we propose is the following 


eq:model 


(43) 


' d t p+ d x {pu) = 0, 

d t pi + d x (p 1 u) = — (p-p{)(p- p 2 )f(p 2 \pi), 

1 £ 

dtP 2 + d x (p 2 u) = -(p - Pi)(p- P 2 )f{pi\p 2 ), 

„ d t (pu) + d x (pu 2 + aip(pi) + a 2 p(p 2 )) = 0, 


where e is a relaxation parameter that determines the rate at which the chemical 
potentials and pressures of the two phases reach equilibrium. The chemical poten¬ 
tial p and the pressure p follow the van der Waals model (11). The source terms 
on the partial densities equations are exactly those of (36) , and involve the relative 
free energy f(pi\pj) which is defined in (34). System (43) is supplemented with 
initial conditions on the velocity u and on the densities p and pi, i = 1,2 satisfying 
the assumption (Hi). 

Combining the mass conservation equation (43)-l and the equations on the par¬ 
tial densities pt, i = 1,2, one can compute the equation satisfied by the volume 
fraction o>\ 


eq:eq_alpha 


(44) 


d t ai + ud x ai 


1 (P~ Pl)(P~ P 2 ) 


[uif(p 2 \pi) - a 2 f(p l \p 2 )\. 


e 


Pi - P2 
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sec:hyperbolicity 


eq:nrj _math 


eq:comp_bal_law 


From (43) one can also recover the two equations on the partial masses atpt and 
deal with a system of the classical form 

' d t p+ d x {pu) = 0 , 


dt(a\P\) + d x (a 1 p 1 u = —— [aip 2 f(p 2 \pi) ~ a 2 pif(pi\pi )\, 

£ Pi — P 2 

dt{a 2 p 2 ) + d x (a 2 p 2 u = — Pl ^ P —— [a 1 p 2 f(p 2 \p 1 ) - a 2 pif(pi\p 2 )}, 

£ Pi ~ P2 

.d t (pu) + d x (pu 2 + aip(pi) + a 2 p(p 2 )) = 0 . 


In the present paper we choose to focus on the system (43), which allows us to 
define the densities pt even when cq = 0 , which is not the case in the last system. 


Remark 4.1. An interesting feature is that the system boils down to the classical 
p— system in pure phases that is when a\a 2 = 0 , including the metastable regions. 


4.2. Hyperbolicity and entropy for the homogeneous system. We introduce 
the mechanical energy 

puf 

(45) £(p,p!,p 2 ,u) = — + F{p,PhPi), 

where the total Helmoltz free energy T is defined by (30). The first result we obtain 
is the decrease in time of this energy. 

Proposition 4.1. The functions, defined in (45), satisfies the following equation 

(46) d t £ + d x (u(£ + aip(pi) + a 2 p{p 2 )) < 0. 

Proof. On the one hand, using the notation (31), one has 

d t F(v) = V r T{v)d t v = -V r F{r)d x (uv) + V r ^(r) • Q, 

where Q = ^ (0, -{p - pi){p - p 2 )f{p 2 \pi), (p - pi)(p - p 2 )f(pi\p 2 )Y- Hence it 
comes 

d t T{r) = -ud x F{ r) - V r J r (r) • r d x u + V r J r (r) • Q. 

Now the expression of V r J r (r) given in (33) leads to 

V r .F(r) • r = ---(((/(pi) - f(p 2 ))p + cxipif(p 2 \pi) - a 2 p 2 f(p 1 \p 2 )). 

Pi - Pi 

Accounting on the definition of the relative free energy (34), it yields 
V r J 7 (r) • r = F(r) + aip(pi) + a 2 p{p 2 ). 


Thus one has 

d t lF(r) = - ud x J r (r ) - (J"(r) + Oip(pi) + a 2 p(p 2 ))d x u + V r J r (r) • Q. 
= ~d x {F{r)u) - (aip(pi) + a 2 p{p 2 ))d x u + V r J r (r) • Q. 

On the other hand, a classical Euler type computation gets 


d 


pu 


+ d x 


pu 


ud x (aip(pi) + a 2 p(p 2 )) = 0 . 


Combining the previous two relations gives 

d t S + d x (u(£ + aip(pi) + a 2 p(p 2 ))) = T{r) Q < 0, 
where the final inequality follows again from the expression (33) of V r J r (r). 


□ 


However, the function £ is not convex everywhere, so that it cannot be considered 
as a mathematical entropy. Indeed £ is convex where T is, and we have the following 
result. 









thm:convex_F 


eq: syst_prim 


eq: AY 


eq:ABC 


eq:speed_c 


eq:eigval 


eq:hyp_cond 
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Theorem 4.1. The total Helmholtz free energy T defined by (30) is convex for 
p £]0, 3[, u € R. and 

• (Pi,P2) e (0,p) x ((p,p“) U (p + , 3)), if p< p~ 

• (pi,P2) G (0,p-) X (p+,3), if p € (p~,p + ) 

• (pi,P2) e ((o ,P~) X ( p + ,p )) x (p, 3), if p> p+ 

Proof. The function T is a convex combination of /(pi) and /(P2) where / is the 
intensive Helmholtz free energy (8). By definition of p_ and p + (see Figure 1), / 
is convex on (0, p~] U [p + , 3), so the result follows. □ 


We turn now to the determination of the eigenvalues of the homogeneous sys¬ 
tem (43). If we set Y = (p, pi, p 2 , u), for smooth solutions, the homogeneous system 
can be written as 


(47) 


d t Y + A(Y)d x Y = 0, 


where the matrix A(Y) is defined by 


(48) 


H(Y) 


( u 

0 

0 

P\ 

0 

u 

0 

Pi 

0 

0 

u 

P2 

W Y ) 

MY) 


u) 


and 


(49) 


Ai(Y) 

M Y) 

4*(Y) 


p(pi) -p(P2) 

p(pl - P 2 ) 

/ Ql - n(p(P 2) -p(pi)) + —p'(pi), 

P(P 1 - P 2 ) p 

-c(p(P2) -p(pi)) + — p'(P2). 

P(Pl - P2) P 


The characteristic equation of A(Y) is given by 


(■u — A) 2 (zi — c — A)(u + c — A), 


with the speed of sound 


(50) c := c(r) = \j ~ («i(P:Pi,P 2 )piP'(pi) + a 2 (p,pi, P 2 )P 2 P 1 (P 2 ))■ 


Thus we obtain three distinct eigenvalues for the matrix Al(Y): 


(51) Ai(Y ) = u-c, A 2 (Y) = A 3 (Y) = u, A 4 (Y) = u + c. 

Note that the eigenvalues are real if r satisfies 

ai(r)pip'(pi) + a 2 (r)p 2 p\p 2 ) > 0. 

Accounting on relations (9), it is equivalent to the following hyperbolicity condition 

(52) ai(r)p?p'(pi) + a 2 (r)pjp’{p 2 ) > 0. 

The right eigenvectors rfi Y), i = 1,..., 4, that satisfy A(Y)r^(Y) = Ai(Y)r*(Y) 
can be chosen as 
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eq:right_eigvect 


sec:struct-waves 


(53) 


n(Y) = 


c 




/-M 


c 

_P 1 


Ai 


Ai 


Pi 

c 

, MY) = 

0 

, MY) = 

1 

, M Y) = 

C 

_P2 


1 


0 


P2 

c 


\ 0 ) 


\ 0 J 


C 


V 1 / V 1 / 


where the quantities Ai, A 2 and A 3 are defined in (49). 

If the densities p, pi and p 2 satisfy (Hi) and (52), the matrix A(Y) is diago- 
nalizable in R and its eigenvectors span the whole space R 4 so that the system is 
hyperbolic. 


4.3. Structure of the waves. In this paragraph we study the structure of the 
waves. Assuming that the densities p, p\ and p 2 satisfy (Hi) and (52), one can 
observe that the waves are either genuinely non linear or linearly degenerate. 

Straightforward computations lead to the following property which will be useful 
in the sequel. 



eq:p_c 



eq:grad_c 


Proposition 4.2. The speed of sound c, function of state r, satisfies the following 
properties relations 

(54)Vp(r) ■ r = pc 2 (r), 


/ 


2c(r) 


c r 


+ 


1 


P P(Pl ~ P2) 

M r ) / P2P'(P2) ~ PlP (Pl) 
P V Pl~P2 

M r ) / P2P'(P2) - pip'(pi) 

V P V Pi - P2 


(PlP'(Pl) - P2P'(P2)) 


+ P'(Pl) +PlP"(Pl) 
+ P'(P2) A P2P”(P2) 


\ 


/ 


Let us start with the waves associated to the wave speed u — c and u + c. 


Proposition 4.3. The characteristic fields associated to the waves speed Ai(Y) = 
u — c and A 4 (Y) = u + c are genuinely non linear i.e. VyAi(Y) • rr(Y) ^ 0 and 
VyA 4 (Y) • 7*4(Y) ^ 0 for admissible state vector Y that is for densities (p,pi,p 2 ) 
satisfying (Hi) and (52). 


eq:gen_non_lin_field 


Proof. We introduce the notation D( r) = (c(r)) 2 . We consider the wave associated 
to the eigenvalue A X (Y). One has 


V Y Ai(Y)-n(Y) = 


(56) 


V r c( r ) 

1 

2c(r) 2 

1 

2pc(r)‘- 


■ri( Y) 

dD dD dD \ 

P ^ +Pl Wi +P2 W2) +1 

(ai(v)p\p"(pi) + a 2 (v)p 2 2 p"(p 2 )) + 1. 


The densities are assumed to be strictly positive. Under the hypothesis (Hi) and 
(52) the mass fractions ai are positive and the second derivative of the van der Waals 
pressure (11) is a strictly negative function of the density. Thus VyAi(Y) •r 1 (Y) / 
0. Similarly we can state that VyA 4 (Y) • 7-4 (Y) ^ 0 that conclude the proof. □ 


We now study the wave associated to the speed u. 



























eq:lin_deg_field 


eq:edp_fraction 


eq:RI_u 


eq:RI_u+c 


eq:RI_u_p 


sec:dom_hyp 


def:dom_inv 
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Proposition 4.4. The characteristic fields associated to the waves A 2 (Y) = Aa(Y) = 
u linearly degenerate i.e. VyA 2 (Y) • r 2 (Y) = 0 and VyAs(Y) • ra(Y) = 0 for ad¬ 
missible state vector Y that is for densities (p, pi,p 2 ) satisfying (52). 

Proof. We deduce from the eigenvalues (51) the relation 

(57) V y A,(Y) • r*(Y) = (0,0,0,if • rfY), 

for i = {1,2}. Then introducing the right eigenvectors (53) in (57), it is easily 
checked that VyA,;(Y) • r*(Y) = 0 for i = {1, 2} and this complete the proof. □ 

We now address the determination of the Riemann invariants of the system. 
These computations are made easier using the following property. 

Proposition 4.5. The mass and volume fractions ol\ and p±, defined by (17) 
and (18), satisfy the following non conservative equations 

d t a\ + ud x ai = 0 , 

^ d t p i + ud x ipi = 0. 

Proposition 4.6. The Riemann invariants associated to the wave of speed u are 

(59) {w,p}, 

with p(r) = ai(r)p(pi) + a 2 (r)p(p 2 ) ■ The volume and mass fractions are Riemann 
invariants associated to the wave of speed u ± c: 

(60) {oti,pi}. 

Proof. Because the field associated to the speed u is linearly degenerate, u is clearly 
a Riemann invariant for this wave. Using the gradient of p with respect to r, a 
straightforward computation gives 

(61) Vyp(Y) • r a (Y) = 0. 

On the other hand the volume fractions on and the mass fractions pt satisfy the 
equations (58). Thus the fractions are Riemann invariants for the waves of speed 
u±c. □ 

The characterization of the third Riemann invariant for the waves of speed u ± c 
is more intricate and is not addressed here. 

4.4. Invariant domains of hyperbolicity for the relaxed system. According 
to Section 4.2, it is clear that the homogeneous system (43) is hyperbolic if and only 
if the densities p, pi and p 2 satisfy (Hi) and the constraint (52) on the speed of 
sound. Nonetheless we are interested in the study and the numerical approximation 
of the whole relaxed system (43), that is taking into account the relaxation term 
with a finite relaxation parameter £ > 0. Actually the domains of hyperbolicity 
of (43) strongly depend on the attraction basins of the dynamical system (36). In 
the present section, we introduce the notion of invariant domains in the same spirit 
as in [9] for diffusive systems. We show that invariant domains U of hyperbolicity 
for the relaxed system (43) are subsets of the attraction basins of the dynamical 
system (36). First note that the hyperbolicity of the homogeneous system (43) solely 
depends on the densities r (t,x) = (p, pi, p 2 ) t (t, x), V(t,x) G R + x R, according to 
the constraint (52) on the speed of sound and not on the velocity u(t, x). Hence we 
consider the following definition of an invariant region. 

Definition 4.1. Let H = {r = (p,pi,p 2 ) e]0, 3[ 3 | 0 < p\ < p < p 2 and p\ < p 2 } a 
subset of the phase space (p,pi,p 2 ) with a Lipschitz continuous boundary dTl. The 
region Q is said to be a invariant domain if 

{fix e R, r(0, x) € H} <t=> Vf > 0, {Vai € R, r (t, x) € fl}. 
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eq:indicS 


eq: SJ 


Stokes 


eq:lcS 


(62) 


We now define some kind of indicator function for such a domain f l: let S be 
defined by 

S : ]0,3[ 3 -4 R 

r=(p,pi,P 2 ) ps{p 1 /p,p 2 /p), 

where s(/3i, /? 2 ) = 1 — l{ / g-< l g 4 < l g+}- Obviously we have 

S(r) = 0 r G tt. 

Next we introduce the nonnegative quantity J 

J : R+ -A R 

1 ^ fR S ( r (t’ X )) dx - 


(63) 


Proposition 4.7. Consider Q a subset of the phase space (p, pi, p 2 ) with a Lipschitz 
continuous boundary dfl and the associated function S defined by (62). Then one 
has the following properties: 

(1) In the sense of distributions we have 

(64) (VS, <j>) = [5] f n(f>(a) da, 

Jd n 

where da is the surface measure on dfl, n the outer normal of Q and 
[S] = S ou t — Si n is the jump of S across the boundary dfl. 

(2) The function S is positively homogeneous of degree 1 so that it verifies the 
Euler relation S(r) = V r S(r) • r. 

(3) The function S satisfies 

(65) (V r S(r), dtX + d x (ur)) = d t S(r) + d x (uS(r)). 


Proof. The first item is a consequence of the Stokes theorem. By construction the 
function S is positively homogeneous of degree 1. Then it satisfies the Euler relation 
given in the second item. Finally following the same steps as in the energy estimate 
(46), we have 

V r S(r)<9 t r + V r S(r)d x (ur) 

= d t S( r) + ud x S{ r) + V r S(r) • r d x u 
= d t S( r) +d x (uS(r)), 

where we use the above Euler relation for S to obtain the last equality. □ 


prop:caract_invdoml 


lem:dom_inv 


We now relate the definition of an invariant domain to the functions S and J 
through several propositions. 

Proposition 4.8. The domain fl is an invariant region if and only if 
{J(0) = 0 => Vf > 0, J(t) = 0}. 

The proof of the Proposition 4.8 relies on the following Lemma. 

Lemma 4.1. LetCi be a subset of the phase space (p, p\, p 2 ) and S defined by (62). 
Then one has 


Vir € R, r(.,a;) G fl J(.) 


5(r(., x))dx = 0. 


The proof of the Lemma relies on the definition of S and its positivity. 

Proof of Proposition j.8. According to the Lemma 4.1 and by the definition of J 
(63) of the quantity J, it follows 

Vx e R, r(0, x) G VL J(0) = 0. 

Using the Lemma 4.1 again, one gets 

Va; G R, r{t, i) G U o J(t) = 0. 
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Combining these two equivalences leads to the conclusion. 


□ 


prop:J_inv 


Proposition 4.9. Let LI be a subset of the phase space (p, pi,p 2 ) and J given by 
(63). Assume J is differentiable. Then it follows 


s J(,) s 0 


LI is an invariant domain. 


Proof. Assume that J(0) = 0. By assumption on the time derivative of J, J{t) < 
J(0) and J(t) > 0 by positivity. Thus J(t) = 0, Vi. Hence according to Proposition 
4.8, the domain LI is invariant. □ 


Corollary 4.1. Let LI be a subset of the plan (p, pi,p 2 ) and S and J the associated 
functions given by (63) and (62). Denote 

Q = 1 (0, ~(p - Pl)(p - P2)f(P2\Pl), (P - Pl)(p ~ P2)f(Pl\P2)f 

£ 

the right-hand side of the relaxed model (43). Then for any r £ LI such that 
lim r(.,x) = lim r(.,x), one has the following assertions 

x —>-+oo x—t — oo 

(Q, Vj-iSYr)) <0 => — Jit) <0 => LI is an invariant domain, 

dt 

Proof. Since Q is the right-hand side of the relaxed model (43), it yields 
(Q, V r S(r)) = (d t v + d x (ur),V r S(v)). 

According to (65) it follows that if (Q, V r S(r)) < 0 then 

(d t r + d x {ur), V r 5(r)) < 0. 

Integrating the above inequality on R gives f R d t S{r)dx < 0, that is 

Proposition 4.9 now leads to the conclusion. □ 


prop:graphf 


Hence in order to check that LI is an invariant domain, one has solely to verify 
that 

<Q,V r S(r))<0. 

Taking the scalar product of VS 1 with the right-hand side Q of the relaxation system 
we obtain 

<Q,V P S(r)) = ^(p-pi)(p - P 2 )( - d Pl Sf{p 2 \pi) + d P2 Sf{p 1 \p 2 )). 

Using equation (64), we obtain formally 

(Q,V r S(r)) = l(p- pi)(p - p 2 )p{ ~ n 1 f{p 2 \pi)+n 2 f{pi\p 2 )), 

where n p = (ni,n 2 ) is the outer normal of the domain Ll p = {(pi, p 2 )\(p, pi, P 2 ) € 
LI}. Note that the domain Ll p are rectangles in the phase space [p\, p 2 )- Now check¬ 
ing the sign of (Q, V r S(r)) is quite straightforward on each part of the boundary 
d Ll p . 

To characterize the invariant regions we study the sign of the relative Helmoltz 
free energy /(.|.) defined by (34). According to Figure 8, one can determine the 
sign of /(.|.) according to the following proposition. 

Proposition 4.10. Let 6 e]0,/4([ and p^ £ [P2 j3 [ such that f{poo\p~) = 0. Then 
the relative Helmoltz free energy satisfies 

• f(p\S) > 0, Vp e]5,Poo[, 

• /(p|p _ ) > 0 (resp. < D), Vp G]5,p“[ (resp. p £]p~,poo[), 
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m 



fig:graphf 


Figure 8. Some reference points on the graph of the Helmoltz 
energy f(p). The blue curve is a sketch of the graph of f{p). The 
green line is the tangent to the graph of / at the point 5. The 
red lines are the tangents of the blue curve at p~ and p + . The 
dashed line is the tangent of the graph of / at p^ defined by 
f(poo,P~) = 0. Depending on the position of the tangent line to 
the blue curve at a given point p , one can determine the sign of 
/(.|p), see Proposition 4.10. 



fig:Apriori 


Figure 9. A priori estimate for (pi,p 2 ): the blue area refers to 
nonattainable states according to (Hi), the green zone is an invari¬ 
ant domain, providing in particular that void cannot appear. 


• f(p\p + ) < 0 (resp. > 0), Vp €]<5,p+[ (resp. p e]p+PooL>, 

• f(p\Poo) > 0, Vp e]<5,poo[. 

Proof. The sign of /(.|a) for any remarkable density a depends on the positition of 
the tangent to the graph of / at the point a, see Figure 8. If the tangent at the point 
a is below the curve (resp. above), one has f(.\a) = /(.) — /(a) — f'(a)(. — a) > 0 
(resp. < 0). □ 

We first state a global a priori estimate which ensures that if (Hi) is satisfied 
at t = 0 then it is preserved for any time by the relaxed system (43). 





















26 


FRANQOIS JAMES AND HELENE MATHIS 


prop:aprioriestimate 


Proposition 4.11. Let 8 £]0,pj] and p^ £ [/?2>3[ (see Figure 8). Then, for any 
0 < p < 3, the domain 


:= {(pi,P2) e {S,p) x ( p,Poo )} 


is invariant. 


Proof. One has to check the sign of (Q, V r S( r)) on each side of the green rectangle 
domain, see Figure 9. On the sides {pi = p} and {p 2 = p} of the rectangle, Q 
vanishes. Now for the side {p\ = <$}, the outer normal is n p = (—1, 0) and 

(Q,V r S( r)) = -{p- 5)(p- p 2 )pf(p2\p~). 

£ 

The product (p — 5)(p — p 2 )p is nonpositive and 

f{P2\p~) = f(p2) - f(P~) - f'(P~)(p2 ~ P~) > 0, 

for p 2 > P 2 , see the green tangent line on Figure 8. A similar argument involving 
/(/’llPoo) works for the side {p 2 = Poo}- □ 


We turn now to determine the invariant domains of hyperbolicity of the relaxed 
model (43) depending on the value of the density p. 


prop:invdom 


Proposition 4.12. Fix 0 < 8 < p\ and let p^ he such that f(poo\p ) = 0. The 
following subsets f l p are invariant domains of hyperbolicity : 

(1) Pure gaseous stable zone. For any 8 < p < p\, 

'■= {(Pi> P 2 ) e]<5,P"[x]p + ,Poo[}, 


see Figure 3. 

(2) Pure liquid stable zone. For any p^ > p > p%, 

’•= {(pi)P2) e]p + ,p[x\p,Poc[}. 

(3) Metastable zones. For any 8 < p < p^, 

^p : = {(pi; P 2 ) £]i5,nhn (p, p~)[x\max(p, p + ), p x [}, 


see Figure 4 for the metastable vapor zone. 

(4) Spinodal zone. For any p~ < p < p + , 

^p : = {(pis P2) €] 5 , p _ [x]p + , poo [}, 


see Figure 2. 


Proof. We only give the proof for the spinodal zone. Following the proof of Propo¬ 
sition 4.11 one has to check that {Q, V r S(r)) on the boundary of the green domain 
Lip of Figure 2. On the side {p\ = p~}, the outer normal is n p = (1,0)*. Thus one 
has 

(Q,V r 5( r)} = l(p-p-)(p-p 2 )p(-f(p 2 \p~)). 

The product (p — p~)(p — P 2 )p is nonpositive since p > 0, p > p~ (because p is 
fixed in [p _ ,p + [) and p 2 > p thanks to hypothesis (Hi). Moreover /(p 2 |p _ ) < 0 for 
any p 2 £]p _ ,Poo[ according to Proposition 4.10. Indeed the red tangent line of the 
blue curve of / at p~ is above the graph of /. Hence (Q, V r S(r)) < 0 on the side 
{pi = p~}. The same kind of arguments are used to prove that (Q,V r S(r)) < 0 
on the three other sides of Ll p . □ 


Note that this characterization of invariant domains of hyperbolicity is formal 
since it relies on the smoothness of the densities r. A possible way to generalize to 
weak solutions is to follow the definition of Hoff [20] . 
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sec:numer-appr 


sec:defin-splitt-strat 


eq:cauchy 


eq:cauchy_disc 


eq:convective 


eq:relax 


eq:W_i 


eq:rusanov 


5. Numerical approximation 

This section is devoted to numerical experiments. We do not wish to elaborate 
here on efficient numerical schemes for this problem, but merely to illustrate some 
typical behaviours of the model. Hence we limit ourselves to a simple finite volume 
scheme, coupled to a time-splitting method for the source terms. Considering the 
stiffness of the problem, a complete numerical study is mandatory but far beyond 
the aim of this paper. We emphasize that we did not implement any specific 
strategy for the non hyperbolicity of the homogeneous system. However for all the 
cases we present, the computed sound velocity is real. This does not prevent from 
possible losses of hyperbolicity, probably due to lack of convergence in the source 
term treatment, see in particular sections 5.2.3 and 5.2.4. 


5.1. Definition of the splitting strategy. We rewrite the system (43) in a more 
compact form, considering the following Cauchy problem 


d t w + d x F{w) = s(w), 

W(t = 0,x) = Wo(a;), Vx £ R, 


where 

W = (p, pi,p 2 ,pw) T , 

F{W) = (pu, piu, p 2 u , pu 2 + aip(pi) + a 2 p(p 2 )) T , 

S(W) = ^(0, -(p - pi){p - p 2 )f(p 2 \pi), (p - pi)(p - p~2)/(pi | p 2 )10), 

and e is the relaxation parameter. Note that we exclude pure phase initial data, so 
that the equations on the partial densities are not multiplied by a* anymore. 

Convective terms and source terms are taken into account by a fractional step 
approach. We denote At the time step and Aa: the length of the cell (cc.j_i/ 2 , x i+ i/ 2 ) 
on the regular mesh. Let W n be the Finite Volume approximation at time t n = 
nAt, n £ N. The approximated solution W n+1 of the Cauchy problem 


(67) 


d t W + d x F{W) = S(W), t £ (t n ,t n+1 ), x £ R, 
W{t n ,x) = W n {x), Vz £ R, 


is approximated by splitting the problem in two steps. The first one corresponds 
to the convective part 


( 68 ) 


d t W + d x F(W) = 0, t £ (t n ,t n+1 ), x £ R 
W{t n ,x) = W n (x), \/x £ R, 


which provides W n ’ . The second steps corresponds to the relaxation process 


(69) 


d t w = S{W ), t £ (f n ,<” +1 ), X £ R 
W{t n ,x) = W n ~(x), Vx £ R, 


which finally gives W n+1 . 

Numerical scheme for the convective part. We consider a classical HLL 
numerical flux. We adopt the following classical notation 

(70) Wp = —*— [ W(t n ,x)dx , n > 0, * £ R. 

Aa; L. 

J x i- 1/2 

The scheme is the following 


(71) 


A x(W?’~ - W.n + At(F? +1/2 - F?_ 1/2 ) = 0, 
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together with 


771/1 

r i+ 1/2 


p(wn, 

srF(W?) - s L F(W ? +1 ) + s L s R (W ? + 1 - W?) 
sr - sl 

F{W? + 1 ), 


i/0 < s L , 

*/sl < 0 < s fi , 
i/0 > sh, 


where s# = max(u" + c7,w" + i + c" + i) and sl = min(u™ - c™,u ™ +1 - c" + i)- The 
time step is subjected to the classical CFL condition 


eq:CFL 


(72) 


At 

Ax 


|A 


max 


< 1, 


where A max is the maximal speed of wave computed on each cell of the mesh. 

Numerical treatment of the source terms. The initial condition for this 
step is W n ~ which is assumed to be admissible, that is {p n ’~,Pi’~,P 2 ~) £ F. The 
total density p and the momentum pu remain unchanged during this step. Only 
the densities p\ and p2 may vary which leads to the following system 


eq:source 


sec:numerical-results 


dtp = d t {pu) = 0, 

(73) dtPi = -~{p - Pi)(p- P2)f(p2\pi), 

dtP -2 = ^{p~ Pi){p~ P 2 )f{pi\p 2 )- 

At this stage we merely use a classic explicit order 4 Runge-Kutta method to 
integrate the source term. 

Such a treatment enforces tough constraints on the time step: the computations 
were performed with 1000 iterations using a time step of 10~ 6 . We emphasize that 
this does not ensure the actual convergence to the equilibrium state. This pleads 
for a more efficient method, for instance a semi-implicit scheme in the spirit of [7]. 

5.2. Numerical results. We present here numerical results that assess the ability 
of the model (43) to capture phase transition and metastable states. 

We consider the van der Waals pressure in its reduced form (13) at a constant 
subcritical temperature T = 0.85. At this fixed temperature the extrema p~ and 
p + of the pressure and the values p\ and p\ defined by the Maxwell construction 
on the chemical potential are given in Table 1. 

We propose test cases with Riemann initial conditions that is 


eq:CI_RP 


(74) 


W(t = 0, x) = Wq ( x ) 


W L , if x < 0, 
W R , if x > 0. 


c:riemann-problem-with 


The following test cases are set on the domain [0,1], with an uniform mesh of 
10000 cells and Neumann boundary conditions at x = 0 and x = 1. 

5.2.1. Riemann problem with phase transition. The initial state Wr and Wr are 


Pl = Pi,l = 0.3, p 2 ,L = P*2 , ul = 0, 
PR = p2,R = 1-9, p\ t R = P*, Ur = 0. 


The left state is a pure stable gas and the right state is a pure stable liquid. Various 
value of £ are considered. The solution is at time t = 0.1s. One can observe the 
appearance of a mixture zone on both sides of the interface, see in particular the 
pressure profile Figure 10-bottom left. The results with e = 10~ 4 or 10 -6 are similar 
expect on the velocity profile (see Figure 10-bottom right) where the intermediate 
state on the left of the interface is slightly modified for e = 10 -4 . 
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fig:transition! 


Figure 10. Riemann problem with phase transition. From top 
left to bottom right: density p, densities rho± and p 2 , volume 
fraction or, speed of sound c, pressure p and velocity u. 


sec:avit-double-raref 


sec:bulle-double-choc 


5.2.2. Cavitation by double rarefaction. The test consists in a liquid state submitted 
to a double rarefaction wave. The initial state is given by p = p\ — 10 -3 , P2 = P21 
pi = p\ and the velocities are ur = 0.3 = — ul • The total density corresponds 
to a metastable liquid, and the initial volume fraction is or ~ 0.000672, which 
means that phase 2 is predominant. The solution is computed at time t = 0.1s. 
We observe on the plot of the volume fraction (Figure 11-second line left) that a 
bubble of stable vapor appears around the interface x = 0.5. The value of e does 
not modify the profile of the bubble. However the pressure profile is sharper for 
e = 10 -6 than for 10~ 4 , see Figure 11-bottom left. 

5.2.3. Nucleation by double shock. The test consists in a pure stable gaseous state 
submitted to a double shock wave. The initial state is given by p = p\ = 0.3, 
P 2 = 1 and the velocities are ur = —0.3 = — ul- The solution is computed at time 
t = 0.4s. Note that phase 2 is not present initially but is fixed in the spinodal 
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fig:cavitl 


Figure 11. Cavitation by double rarefaction. From top left to 
bottom right: density p , densities rho\ and P 2 , volume fraction 
ot\, speed of sound c, pressure p and velocity u. 


zone. We observe two different behaviours of the solution depending on the value 
of e. For e = 10 -6 , the profile of the volume fraction (see Figure 12-second line 
left) shows that a droplet of liquid appears around x = 0.5 which is not the case for 
e = 10 -4 where there is no droplet. For e = 10~ 6 , one observes on the density plot 
(see Figure 12-top left) that the liquid state inside the droplet admits a density 
p close to P 2 with small oscillations. The droplet is surrounded by two mixture 
areas with p = p\. The pressure curve inside the droplet presents oscillations (see 
Figure 12-bottom left) which might be due to a loss of hyperbolicity due to the lack 
of accuracy in the approximation of the source (63). 


sec:perturb-meta 


5.2.4. Acoustic perturbation of a metastable state. This example consists in a con¬ 
stant metastable vapor state, with a perturbation in the velocity. The initial state 
is p = 0.42, pi = 0.32, P 2 = 0.52 and the velocities are ul = 0.4, ur = 0. Both 
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Figure 12. Nucleation by double shock. From top left to bottom 
right: density p 1 densities rho\ and p 2 , volume fraction aq, speed 
of sound c, pressure p and velocity u. 


densities p\ and p 2 are in the metastable state, and we impose a compression from 
the left with velocity 0.4. 

The compression induces the appearance of droplet of pure liquid which moves 
from the left to the right, see the time evolution on Figure 13. With a smaller 
velocity perturbation the structure of the waves is similar, but there is no creation 
of a droplet at the interface. One can check on Figure 13- top left that the density 
p inside the droplet is larger than p \. The droplet is surrounded by two areas with 
a mixture state with p = p\. The velocity and pressure profiles exhibit spikes on 
both sides of the droplet, which amplitude decreases when e decreases, see Figures 
13, 14 and 15. Notice also that the pressure in the mixture zone is not at the value 
p*. It seems that when e decreases the value is closer, this may indicate that the 
source term has not reached the equilibrium state yet. 
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Figure 13. Perturbation of a metastable state at t=0.1s. From 
top left to bottom right: density p, densities rhoi and p 2 , volume 
fraction aq, speed of sound c, pressure p and velocity u. 


sec:conclusion 


6. Conclusion 

The core of this work is the formalization in terms of a dynamical system of the 
thermodynamics of phase transition, using the van der Waals EoS. It leads to a 
mathematical characterization of metastable states, compared to stable coexistent 
two-phase states. The dynamical behaviour of the solutions is crucial here, and 
seems to preclude the use of instantaneous exchange kinetic. When coupled to a 
simple hydrodynamic model, namely the isothermal Euler equations, it evidences 
abilities to cope with metastable states as well as bubble or droplet generation. 
This preliminary study gives rise to a wide bunch of open questions and problems, 
and we believe that the methodology can be used in a much larger context. 

First, in the same isothermal context, the construction of the dynamical system 
(the right-hand side in the extended Euler equations) can be addressed. We de¬ 
liberately used a simple and readable function, which possibly could be improved. 
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Figure 14. Perturbation of a metastable state at t=0.2s. From 
top left to bottom right: density p, densities rhoi and p 2 , volume 
fraction a 1; speed of sound c, pressure p and velocity u. 


In any case, the behaviour of the coexistence zone around the interface has to be 
investigated in more details, as well as the role of e. 

Next, an obvious mandatory issue is the numerical treatment of the coupled 
system. We have chosen here the simplest numerical strategy that allowed us to 
illustrate our purpose. The explicit treatment of the stiff relaxation term enforces 
tough constraints on the time step, and prevents the simulation of more realistic 
metastable cases. 

Finally, we attend to include temperature dependance to obtain a fully heat, 
mass and mechanical transfer model in order to compare our results to those of [32] 
and [35]. 


knowles91 


References 

[1] R. Abeyaratne and J. K. Knowles. Kinetic relations and the propagation of phase boundaries 
in solids. Arch. Rational Mech. Anal., 114(2):119—154, 1991. 
































































fig:meta4 


faccanoni07 


chalons09 


BN86 


BH05 


bedjaoui05 


IMEX09 


callen85 


34 


FRANQOIS JAMES AND HELENE MATHIS 





X 




X 



Figure 15. Perturbation of a metastable state at t=0.4s. From 
top left to bottom right: density p, densities rho\ and p 2 , volume 
fraction aq, speed of sound c, pressure p and velocity u. 
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